#include<stdio.h>
#include<Arcroll.h>
#include<Algorithm/points_cloud_algorithm.h>

#include<tool/libcell_tools_view.h>

#include<libcell_macros.h>
#include<tool/libcell_tools_algorithm.h>

#include<Algorithm/subdivision_of_polygon.h>
#include<Algorithm/threeD_elementary_geometry_tool_library.h>

//#include<tools/tools_heap.h>

#include<time.h>


#define quote lib_cell_quote

static inline double compute_visual_avarage_distance(Node* n)
{
	double* box=mesh_compute_bounding_box_from_node(n);
	double min1=box[3]-box[0],min2=box[4]-box[1],min3=box[5]-box[2];
	double mid=min1>=min2?(min3>=min1?(min1):(min2>=min3?(min2):(min3))):(min3>=min2?(min2):(min1>=min3?(min1):(min3)));
	double vdis=mid/70.0;
	return vdis;
}
typedef struct Interactive_Show_Info{
	Node* cur_n;
	KD_Node*kdtree;
	double r;
	Node*marked_vs;
	Mesh*m;	
	Viewer_Something* vs;
	Viewer_Opengl_Interpreter* voi;

}Interactive_Show_Info;
static inline void interactive_show_info_init(Interactive_Show_Info* isi)
{
	isi->cur_n=NULL;
	isi->kdtree=NULL;
	isi->r=0;
	isi->marked_vs=NULL;
	isi->m=NULL;
	isi->vs=NULL;
	isi->voi=NULL;

}

static inline void interactive_show_info_free(Interactive_Show_Info* isi)
{
	if(isi==NULL)
	{
		return ;
	}
	free_node(isi->marked_vs);
	free(isi);
}

void interactive_show_info_key_callback(Viewer_Intera*vi)
{
	Interactor_GlobalInfo* g_info=vi->g_info;
	Interactive_Show_Info* isi=(Interactive_Show_Info*)(vi->representation);
	if(isi->cur_n==NULL)
	{
		return ;
	}
	int mark=0;
    if(g_info->key_action==VIEWER_PRESS&&g_info->key==VIEWER_KEY_T)
    {
    	isi->cur_n=(Node*)(isi->cur_n->Next);
    	if(isi->cur_n==NULL)
    	{
    		return ;
    	}
    	for(Node* nit=isi->marked_vs;nit!=NULL;nit=(Node*)(nit->Next))
		{
			template_v* v=(template_v*)(nit->value);
			v->prop=NULL;
		}
		free_node(isi->marked_vs);isi->marked_vs=NULL;
		template_v* v=(template_v*)(isi->cur_n->value);
		isi->marked_vs=kd_tree_find_nearest_sphere(v,isi->r,isi->kdtree);
		for(Node* nit=isi->marked_vs;nit!=NULL;nit=(Node*)(nit->Next))
		{
			template_v* v=(template_v*)(nit->value);
			v->prop=&mark;
		}	
		int i=0;
		Viewer_Points*vp=(Viewer_Points*)(isi->vs->evolution);
		for(auto vit=isi->m->v_begin(isi->m);vit!=isi->m->v_end(isi->m);vit++)
		{
			vp->color[i*4+0]=0.0;vp->color[i*4+1]=0.8;vp->color[i*4+2]=0.0;vp->color[i*4+3]=1.0;
   			if(quote(vit)->prop!=NULL)
   			{
   				vp->color[i*4+0]=1.0;
   			}
			i++;
		}
		isi->voi->update_data(isi->voi);
    }
}




Viewer_Something* test_show_mesh_points(Viewer_World* vw,Mesh* m)
{
    Node* n=vw->create_something(vw,"Points");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
    vp->Data_rows=m->num_v(m);
   	vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
   	int i=0;
   	vp->color=(float*)malloc(sizeof(float)*4*vp->Data_rows);
   	//memset(vp->color,1,sizeof(float)*4*vp->Data_rows);	
   	for(auto vit=m->v_begin(m);vit!=m->v_end(m);vit++)
   	{
   		for(int j=0;j<3;j++)
   		{
   			vp->Data[i*3+j]=quote(vit)->point[j];
   		}	
   		vp->color[i*4+0]=0.0;vp->color[i*4+1]=0.8;vp->color[i*4+2]=0.0;vp->color[i*4+3]=1.0;
   		if(quote(vit)->prop!=NULL)
   		{
   			vp->color[i*4+0]=1.0;
   		}
   		i++;
   	} 
    //vp->set_color(vp,0.9,0.5,0.5,1.0); 
    vp->pointsize=10.0;
    free_node(n);
    return vs;
}

void inline test_show_mesh_cells(Viewer_World* vw,Mesh*m)
{
    Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=m->num_v(m);
    vf->Data_index_rows=m->num_c(m);
    get_data_from_2dim_cell(m,&(vf->Data),&(vf->Data_index));
     
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.9,0.6,0.7,1.0);

    vf->normal_rows=vf->Data_index_rows;
    vf->compute_normal(vf);
    
    free_node(n); 

}
static inline void test_show_mesh_lines(Viewer_World* vw,Mesh* m)
{
    Node* n=vw->create_something(vw,"Edges");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    ve->Data_rows=m->num_v(m);
    ve->Data_index_rows=m->num_f(m);
    get_lines_data_from_2dim_cell(m,&(ve->Data),&(ve->Data_index));
    ve->color_rows=ve->Data_index_rows;
    ve->set_color(ve,1.0,0.0,0.6,1.0); 
    free_node(n);
}
static inline void swap(double* values,int left,int right)
{

    double temp=values[left];
    values[left]=values[right];
    values[right]=temp;

}

static void sort_double_value(double* values,int left,int right)
{
    if(left>=right)
    {
        return;
    }
    int m; 
    m = left;  
    for (int i=left+1;i<=right;i++)
    {
        if(values[i]<values[left])  
        { 
            swap(values,++m,i);
        }  
    }   
    swap(values,left,m);  
    sort_double_value(values,left,m-1);  
    sort_double_value(values,m+1,right);
}

typedef struct  Normal_Boundary_Info{
	double normal[3];
	int is_boundary;
	int is_correct_normal;
}Normal_Boundary_Info;


static inline void normal_boundary_info_init(Normal_Boundary_Info* nbi)
{
	memset(nbi->normal,0,sizeof(double)*3);
	nbi->is_boundary=0;
	nbi->is_correct_normal=0;
}
static inline void normal_boundary_info_free(Normal_Boundary_Info*nbi)
{
	if(nbi==NULL)
	{
		return;
	}
	free(nbi);
}

Node* points_cloud_compute_normal_and_get_boundary_points(Node* n,KD_Node*kdtree,int num)
{
	Node* re=NULL;
	MatrixXd A(3,3);
    clock_t t = clock();
 	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
    	Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)malloc(sizeof(Normal_Boundary_Info));
    	normal_boundary_info_init(nbi);
    	template_v* v=(template_v*)(nit->value);
    	v->prop=nbi;
    	Node* n1=kd_tree_find_n_adjacent_verticesn(v,kdtree,num);
    	double* axis=compute_pca_result_using_eigen(n1);	
    	nbi->normal[0]=axis[6];nbi->normal[1]=axis[7];nbi->normal[2]=axis[8];
    	memmove(A.data(),axis,sizeof(double)*9);

    	//printf("%lf %lf \n",nbi->normal[0],A(0,2));
    		
    	MatrixXd B=A.inverse();
    	VectorXd C(3);	
    	double *thetas=(double*)malloc(sizeof(double)*(num-1));
    	int i=0;
    	for(Node* nit1=n1;nit1!=NULL;nit1=(Node*)(nit1->Next))
    	{
    		template_v* v1=(template_v*)(nit1->value);
    		if(v1==v)
    		{
    			continue;
    		}	
    		for(int j=0;j<v1->point_size;j++)
     		{
     			C.coeffRef(j,0)=v1->point[j]-v->point[j];
     		}
     		MatrixXd temp_m= (B*C);
     		double cosx=temp_m(0),sinx=temp_m(1);
     		double dis=SAFE_SQRT(cosx*cosx+sinx*sinx);
     		cosx/=dis;sinx/=dis;
     		double alpha=acos(cosx);
     		
     		thetas[i]=(sinx<0 ? (2*M_PI-alpha):alpha);		
     		i++;
    	}
		sort_double_value(thetas,0,num-2);
		//printf("M_PI:%lf\n",M_PI);
		for(int j=0;j<num-1;j++)
		{	
			double  alpha=thetas[(j+1)%(num-1)]-thetas[j];
			j==num-2?(alpha=2*M_PI+alpha):(alpha=alpha);	
			if(alpha>M_PI)
			{
				nbi->is_boundary=1;
				re=node_overlying(re,v);
				break;
			}
		}
		free(thetas);
    	free(axis);
		free_node(n1);
    }
    printf("运行时间%.5f秒:%d\n", (double)(clock() - t) / CLOCKS_PER_SEC,0); 
	return re;
} 


typedef struct Points_Open_Set{
	//n :: set of vertices
	//adjacent:: bigcap of open sets 
	Node* n,*adjacent;
	template_v* v;
	double * normal;
}Points_Open_Set;

static inline void points_open_set_init(Points_Open_Set* pos)
{
	pos->n=NULL;pos->v=NULL;pos->adjacent=NULL;
	pos->normal=NULL;	
}

static inline void points_open_set_free(Points_Open_Set* pos)
{
	if(pos==NULL)
	{
		return ;
	}
	free_node(pos->adjacent);
	free_node(pos->n);
	if(pos->normal!=NULL)
	{
		free(pos->normal);
	}
	free(pos);
}

// @ n:: is vertices
//@ tree:: return map v to open_sets
// 

Node* compute_points_open_sets(Node* n,Int_RB_Tree* tree,double r,KD_Node* kdtree)
{
	Node* re=NULL;
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v=(template_v*)(nit->value);
		if(tree->find(tree,v->id)==NULL)
		{
			Points_Open_Set* pos=(Points_Open_Set*)malloc(sizeof(Points_Open_Set));
			points_open_set_init(pos);
			pos->v=v;
			pos->n=kd_tree_find_nearest_sphere(v,1.5*r,kdtree);
			re=node_overlying(re,pos);
			for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
			{
				template_v* v1=(template_v*)(nit1->value);
				Node* value=(Node*)(tree->find(tree,v1->id));
				if(value==NULL)
				{
					value=node_overlying(value,pos);
					tree->insert(tree,v1->id,value);
				}
				else{
					Node* values_n=node_reverse(value);
					node_pushback(values_n, pos);
				}
			}			
		}	
		else
		{
			continue;
		}
	}
	printf(" size1 :%d size2:%d\n",tree->size,node_size(n));
	template_v**values=(template_v**)malloc(sizeof(template_v*)*node_size(re));
	int i=0;
	for(Node* nit=re;nit!=NULL;nit=(Node*)(nit->Next))
	{
		values[i]=((Points_Open_Set*)(nit->value))->v;
		i++;
	}	
	KD_Node* kdtree1= create_kd_tree(values,i,0);	
	for(Node* nit=re;nit!=NULL;nit=(Node*)(nit->Next))
	{
		Points_Open_Set* pos=(Points_Open_Set*)(nit->value);
		Node* n1=kd_tree_find_nearest_sphere(pos->v,3.0*r,kdtree1);

		Node* n2=NULL;
		for(Node* nit1=n1;nit1!=NULL;nit1=(Node*)(nit1->Next))
		{
			template_v* v=(template_v*)(nit1->value);
			Node* n3=(Node*)(tree->find(tree,v->id));
			if(n3==NULL||node_size(n3)!=1)
			{
				printf("cuwu\n");
				return re;
			}	
			n2=node_overlying(n2,n3->value);
		}
		pos->adjacent=n2;
		free_node(n1);
	}	
	free_kdnode(kdtree1);
	free(values);
	return re;
}

int  adjust_normal_of_one_points_open_set(Points_Open_Set* pos,double* normal)
{
	if(pos->normal!=NULL)
	{
		return 0;
	}
	pos->normal=(double*)malloc(sizeof(double)*3);
	memset(pos->normal,0,sizeof(double)*3);
	for(Node* nit=pos->n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v=(template_v*)(nit->value);
		Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)(v->prop);
		if(nbi->is_correct_normal==0&& inner_product(nbi->normal,normal,3)<0)
		{	
			nbi->normal[0]=-nbi->normal[0];nbi->normal[1]=-nbi->normal[1];nbi->normal[2]=-nbi->normal[2];	
		}	
		nbi->is_correct_normal=1;
		pos->normal[0]+=nbi->normal[0];pos->normal[1]+=nbi->normal[1];pos->normal[2]+=nbi->normal[2];
	}
	
	return 1;
}
static inline double* tools_compute_normal(template_v*v1,template_v* v2,double* dir)
{

	double dir1[3]={v2->point[0]-v1->point[0],v2->point[1]-v1->point[1],v2->point[2]-v1->point[2]};
	if(normalize(dir1,3)==0)
	{
		printf("cuowu\n");		
		return NULL;
	}
	double* re=(double*)malloc(sizeof(double)*3);
	double fac=inner_product(dir,dir1,3);
	re[0]=dir[0]-fac*dir1[0];re[1]=dir[1]-fac*dir1[1];re[2]=dir[2]-fac*dir1[2];				

	return re;
}

void set_normal_vector_of_manifold_points(Node* n)
{
	if(n==NULL)
	{
		return ;
	}
	double * normal=NULL,*new_normal=NULL;
	Node* n1=NULL;
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		Points_Open_Set* pos =(Points_Open_Set*)(nit->value);
		if(pos->normal==NULL)
		{
			n1=node_overlying(n1,pos);
			new_normal=(double*)malloc(sizeof(double)*3);
			n1->traits=new_normal;
			Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)(pos->v->prop);
			new_normal[0]=nbi->normal[0];new_normal[1]=nbi->normal[1];new_normal[2]=nbi->normal[2];
			break;
		}
	}
	if(n1==NULL)
	{
		return ;
	}
	printf("hrere\n");
	while(n1!=NULL)
	{		
		Node* temp_n=n1;n1=NULL;
		for(Node* nit=temp_n;nit!=NULL;nit=(Node*)(nit->Next))
		{
			Points_Open_Set* pos=(Points_Open_Set*)(nit->value);	
			normal=(double*)(nit->traits);
			int flag=adjust_normal_of_one_points_open_set(pos,normal);
			if(flag==0)
			{
				
			}		
			else
			{
				for(Node* nit1=pos->adjacent;nit1!=NULL;nit1=(Node*)(nit1->Next))
				{
					Points_Open_Set* pos1=(Points_Open_Set*)(nit1->value);
					if(pos1->normal==NULL)
					{
						n1=node_overlying(n1,pos1);
	
						new_normal=tools_compute_normal(pos1->v,pos->v,pos->normal);
						n1->traits=new_normal;
					}	
				}		
			}
			free(normal);	
			nit->traits=NULL;	
		}
		free_node(temp_n);
		printf("size:%d\n",node_size(n1));
	}

	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		Points_Open_Set* pos=(Points_Open_Set*)(nit->value);
		if(pos->normal==NULL)
		{
			continue;
		}	
		for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
		{
			template_v* v=(template_v*)(nit1->value);
			Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)(v->prop);
			if(inner_product(nbi->normal,pos->normal,3)<0)
			{
				nbi->normal[0]=-nbi->normal[0];	nbi->normal[1]=-nbi->normal[1];	nbi->normal[2]=-nbi->normal[2];	
			}		
		}	
	}	
	// 
	//set_normal_vector_of_manifold_points(n);
} 
//  

//
int tools_boundary_lines_is_correct_direction(template_v* v1,template_v* v2,KD_Node* kdtree,double r)
{
	template_v* v=(template_v*)malloc(sizeof(template_v));
	Vertex_init_(v);	
	v->point_size=3;
	v->point=(double*)malloc(sizeof(double)*3);
	v->point[0]=(v1->point[0]+v2->point[0])/2.0;v->point[1]=(v1->point[1]+v2->point[1])/2.0;v->point[2]=(v1->point[2]+v2->point[2])/2.0;
	Node* n=kd_tree_find_nearest_sphere(v,r,kdtree);
	double dir2[3]={ v1->point[0]-v2->point[0], v1->point[1]-v2->point[1],v1->point[2]-v2->point[2] };	
	double normal[3]={0},normal1[3]={0};
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v3=(template_v*)(nit->value);
		Normal_Boundary_Info*nbi=(Normal_Boundary_Info*)(v3->prop);	
		normal[0]+=nbi->normal[0]; normal[1]+=nbi->normal[1];normal[2]+=nbi->normal[2]; 	 		
		if(v3==v1||v3==v2)
		{continue;}
		double dir1[3]={v3->point[0]-v2->point[0], v3->point[1]-v2->point[1],v3->point[2]-v2->point[2]};

		double*dir=out_product(dir1,dir2);
		normal1[0]+=dir[0]; normal1[1]+=dir[1];normal1[2]+=dir[2]; 	 		
		free(dir);
	}
	free_node(n);
	free_Vertex(v);
	if(inner_product(normal1,normal,3)>0)
	{
		return 1;
	}else
	{
		return 0;
	}
	return 1;
}
double* tools_compute_adjacent_normal(Node* n,double r)
{
	double distance=0;
	Node* nit=n;
	template_v* v1=(template_v*)(nit->value);
	nit=(Node*)(nit->Next);
	double * re=(double *)malloc(sizeof(double)*3);
	memset(re,0,sizeof(double)*3);

	while(nit!=NULL&&distance<1.5*r)
	{

		template_v* v2=(template_v*)(nit->value);
		double temp_dis=distance_of_two_points(v1->point,v2->point,3);
		double dir[3]={v1->point[0]-v2->point[0],v1->point[1]-v2->point[1],v1->point[2]-v2->point[2]};
		re[0]+=dir[0];re[1]+=dir[1];re[2]+=dir[2];

		distance+=temp_dis;	
		v1=v2;
		nit=(Node*)(nit->Next);
	}	
	return re;

}

int  increasing_boundaries_points_to_loop(Node**n,Int_RB_Tree*tree,double r ,KD_Node* kdtree1)
{
	Node* node_reversn=node_reverse(*n);
	template_v* reverse_v=(template_v*)(node_reversn->value);
	template_v* v=(template_v*)((*n)->value);
	while(v!=reverse_v)
	{	
		Node* n1=kd_tree_find_nearest_sphere(v,2.5*r,kdtree1);
		if(n1==NULL)
		{
			break;
		}
		//这个法向信息应该是某种邻域内的信息
		double* dir1=tools_compute_adjacent_normal(*n,r);
		normalize(dir1,3);
		template_v* v3=NULL;
		double scores=0; 		
		//double scores=distance_of_two_points(v->point,value[0]->point,3)/3.0*r;
		for(Node* nit=n1;nit!=NULL;nit=(Node*)(nit->Next))
		{
			template_v* v2=(template_v*)(nit->value);
			double dir2[3]={v2->point[0]-v->point[0],v2->point[1]-v->point[1],v2->point[2]-v->point[2]};
			normalize(dir2,3);
			if(v2==v||tree->find(tree,v2->id)!=NULL||(v2!=reverse_v&&node_find(*n,v2)!=NULL))
			{
				continue;
			}
			if(inner_product(dir1,dir2,3)<-0.4)
			{
				continue;
			}
			if(v3==NULL)
			{
				v3=v2;
				scores=inner_product(dir1,dir2,3)+1-distance_of_two_points(v2->point,v->point,3)/(2.5*r);
			}
			else{
				double temp_scores=inner_product(dir1,dir2,3)+1-distance_of_two_points(v2->point,v->point,3)/(2.5*r);
				if(temp_scores>scores)
				{
					scores=temp_scores;
					v3=v2;
				}
			}	
		}
		free(dir1);
		free_node(n1);
		if(v3==NULL)
		{
			break;
		}
	
		
		v=v3;
		*n=node_overlying(*n,v);
	}	
	if(v!=reverse_v||node_size(*n)<=2)
	{
		return 0;
	}
	return 1;	
}
//@kdtree1:: kd tree of boundaries
//  
void connect_boundary_points(Mesh*mesh,KD_Node* kdtree,double r,Node* boundaries, KD_Node* kdtree1,Int_RB_Tree* tree)
{
	printf("boundaries :%d\n",node_size(boundaries));

	int i=0;
	for(Node* nit=boundaries;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v=(template_v*)(nit->value);

		if(tree->find(tree,v->id)==NULL)
		{
			KD_Node* kdnode=kd_tree_find_nearest_and_except_vs(v,kdtree1,&v,1);
			if(distance_of_two_points(kdnode->value->point,v->point,3)>2.5*r)
			{
				continue;
			}
			Node* n=NULL;
			if(tools_boundary_lines_is_correct_direction(v,kdnode->value,kdtree,r)==0)
			{
				n=node_overlying(n,kdnode->value);
				n=node_overlying(n,v);
			}
			else
			{
				n=node_overlying(n,v);
				n=node_overlying(n,kdnode->value);
			}
			//printf("libo\n");
			if(increasing_boundaries_points_to_loop(&n,tree,r,kdtree1)==0)
			{
				free_node(n);
				continue;
			}
			else{
				Node* nit1=n;
				while(nit1->Next!=NULL)
				{
					template_v* v1=(template_v*)(nit1->value);
					template_v* v2=(template_v*) (((Node*)(nit1->Next))->value);
					template_v* vv1=(template_v*)(tree->find(tree,v1->id));
					if(vv1==NULL)
					{
						vv1=mesh->create_vertexv(mesh,v1->point,3);
					}
					template_v* vv2=(template_v*)(tree->find(tree,v2->id));
					if(vv2==NULL)
					{
						vv2=mesh->create_vertexv(mesh,v2->point,3);
					}
					vv1->prop=v1;
					vv2->prop=v2;
					//printf("create cell\n");
					template_v*vs[2]={vv2,vv1};
					mesh->create_cellv(mesh,vs,2);
					tree->insert(tree,v1->id,vv1);tree->insert(tree,v2->id,vv2);
					nit1=(Node*)(nit1->Next);

				}	
				free_node(n);
				i++;
				// if(i>5)
				// {
				// 	break;
				// }
			}	
		}
	}
	printf("mesh num c:%d\n",mesh->num_c(mesh));

}
Mesh* test_connect_boundaries(Node* boundaries,KD_Node*kdtree,double r)
{
	Mesh* mesh=(Mesh*)malloc(sizeof(Mesh));
	Mesh_init(mesh);
	mesh->dimension=1;
	mesh->simplex=1;
	mesh->manifold_require=1;

	Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree);
	template_v ** values=(template_v**)malloc(sizeof(template_v*)*node_size(boundaries));
	int i=0;
	for(Node* nit=boundaries;nit!=NULL;nit=(Node*)(nit->Next))
	{
		values[i]=(template_v*)(nit->value);
		i++;
	}
	KD_Node*kdtree1= create_kd_tree(values,i,0);
	printf("%lf \n",r);
	connect_boundary_points(mesh,kdtree,r,boundaries,kdtree1,tree);
	//printf("mesh num c:%d\n",mesh.num_c(&mesh));
	free_kdnode(kdtree1);
	free(values);
	int_rb_tree_free(tree);
	return mesh;
	//Mesh_free(&mesh);
}
static inline void filled_normal_info_to_new_vertex(template_v* v,double r ,KD_Node*kdtree)
{
	
	//Node* n=kd_tree_find_nearest_sphere(v,2.0*r,kdtree);	
    Node* n=kd_tree_find_n_adjacent_verticesn(v,kdtree,8);
    double *pcas= compute_pca_result_using_eigen(n);
	Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)malloc(sizeof(Normal_Boundary_Info));
    normal_boundary_info_init(nbi);
    v->prop=nbi;
	nbi->normal[0]=pcas[6]; nbi->normal[1]=pcas[7]; nbi->normal[2]=pcas[8];

	// double normal[3]={0,0,0};
	// for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	// {
	// 	template_v* v1=(template_v*)(nit->value);
	// 	if(v!=v1&&v1->prop!=NULL)
	// 	{
	// 		Normal_Boundary_Info*nbi1=(Normal_Boundary_Info*)(v1->prop);
	// 		normal[0]+=nbi1->normal[0];normal[1]+=nbi1->normal[1];normal[2]+=nbi1->normal[2];
	// 	}
		
	// }
	// //normalize(normal,3);
	// if(inner_product(normal,nbi->normal,3)<0)
	// {
	// 	nbi->normal[0]=-nbi->normal[0];	nbi->normal[1]=-nbi->normal[1];	nbi->normal[2]=-nbi->normal[2];	
	// }
	// nbi->normal[0]+=normal[0];nbi->normal[1]+=normal[1];nbi->normal[2]+=normal[2];
	// normalize(nbi->normal,3);
    free(pcas);
    free_node(n);
}
// n is a boundaries
//kdtree is to update 

Node* points_cloud_filled_pre_handle(Mesh*points ,Node* n,double r,KD_Node*kdtree)
{

	Node* nit=n;
	//对于长边进行插值
	while(nit!=NULL)
	{
		Node* nit1=(Node*)(nit->Next);
		nit1==NULL? (nit1=n):(nit1=nit1);
		template_v* v1=(template_v*)(nit->value),*v2=(template_v*)(nit1->value);
		double dis=distance_of_two_points(v1->point,v2->point,3);
		if(dis>r)
		{
			double p[3]={(v1->point[0]+v2->point[0])/2.0, (v1->point[1]+v2->point[1])/2.0, (v1->point[2]+v2->point[2])/2.0};
			template_v*nv= points->create_vertexv(points,p,3);

			printf("one insert %d\n",nv->id); 
			node_insertback(nit,nv);
			kd_tree_insert_vertex(kdtree,nv);	
			filled_normal_info_to_new_vertex(nv,r,kdtree);
			continue;	
		}
		nit=(Node*)(nit->Next);
	}



	Node* re=n;
	return re;
}
void points_cloud_filled_pre_handle1(Mesh*points ,Node** n,double r,KD_Node*kdtree)
{
	template_v** boundaries=(template_v**)malloc(sizeof(template_v*)*node_size(*n));
	int i=0;
	for(Node* nit=*n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		boundaries[i]=(template_v*)(nit->value);
		i++;
	}
	Node* re=NULL;
	re=node_overlying(re,boundaries[0]);	
	double sum=0;
	for(int j=1;j<=i;j++)
	{
		double dis=distance_of_two_points(boundaries[j%i]->point,boundaries[j-1]->point,3);
		sum+=dis;
		if(sum>r)
		{
			re=node_overlying(re,boundaries[j-1]);
			sum=dis;	
		}
	}	

	free(boundaries);

	free_node(*n);
	*n=re;
}
 

//
//
Node*  points_cloud_filled_get_one_boundaries(Mesh* points,Mesh* edges,Int_RB_Tree* tree)
{
	Node* re=NULL;
	template_c*c=NULL;
	for(auto cit=edges->c_begin(edges);cit!=edges->c_end(edges);cit++)
	{
		if(tree->find(tree,quote(cit)->id)!=NULL)
		{
			continue;
		}
		c=quote(cit);	
	}
	if(c==NULL)
	{
		return NULL;
	}
	template_c* cit=c;

	while(cit!=NULL)
	{
		tree->insert(tree,cit->id,cit);
		re=node_overlying(re,(template_v*)(cit->vertices[1]->prop));
		for(Node* nit=(cit->vertices[0]->cells);nit!=NULL;nit=(Node*)(nit->Next))
		{
			if(((template_c*)(nit->value))!=cit)
			{
				cit=(template_c*)(nit->value);
				break;
			}	
		}	
		if(cit==c)
		{
			break;
		}

	}
	return re;
}
Node* points_cloud_filled_intersections_of_boundary(Node*n,double r,double* normal )
{
	if(n==NULL)
	{
		return NULL;
	}
	Int_RB_Tree*tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree)), *tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree);int_rb_tree_init(tree1);
	Int_RB_Tree* tree3=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree3);
	template_v**values=(template_v**)malloc(sizeof(template_v*)*node_size(n)),**temp_values=(template_v**)malloc(sizeof(template_v*)*node_size(n));

	int i=0;
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{

		values[i]=(template_v*)(nit->value);
		//printf("values i:%d %d\n",i,values[i]->id);
		int *value=(int*)malloc(sizeof(int));
		*value=i;
		double* p=(double*)malloc(sizeof(double)*3);
		p[0]=values[i]->point[0];p[1]=values[i]->point[1];p[2]=values[i]->point[2];
		double temp=inner_product(p,normal,3);
		p[0]=p[0]-temp*normal[0];p[1]=p[1]-temp*normal[1];p[2]=p[2]-temp*normal[2];
		tree1->insert(tree1,values[i]->id,p);
		tree->insert(tree,values[i]->id,value);
		i++;
	}
	memmove(temp_values,values,sizeof(template_v*)*node_size(n));
	KD_Node* kdtree=create_kd_tree(values,i,0);

	free(values);
	values=temp_values;

	Node* re=NULL;
	for(int j=0;j<i;j++)
	{
		template_v* vv1=values[j];
		template_v* vv2=values[(j+1)%i];
		Node* n1=kd_tree_find_nearest_sphere(vv1,2.0*r,kdtree);	
		double* p1=(double*)(tree1->find(tree1,vv1->id));
		double* p2=(double*)(tree1->find(tree1,vv2->id));
		
		if(tree3->find(tree3,vv1->id)!=NULL)
		{
			continue;
		}
		for(Node* nit=n1;nit!=NULL;nit=(Node*)(nit->Next))
		{
			template_v* v1=(template_v*)(nit->value);
			if(v1==vv1||v1==vv2||v1==values[(j-1+i)%i])
			{
				continue;
			}
			double* p3=(double*)(tree1->find(tree1,v1->id));
			int k= *((int*)(tree->find(tree,v1->id)));	
			template_v* v2=values[(k+1)%i];
			double* p4=(double*)(tree1->find(tree1,v2->id));

			if(is_line_segment_intersecting_with_line_segment(p1,p2,p3,p4)==1)
			{	
				tree3->insert(tree3,v1->id,v1);
				printf("intersetion vv1 id:%d  %lf %lf %lfv1:%d\n",vv1->id,vv1->point[0], vv1->point[1],vv1->point[2],v1->id);

				template_v**vs=(template_v**)malloc(sizeof(template_v*)*4);
				vs[0]=v1;vs[1]=v2;vs[2]=vv1;vs[3]=vv2;
				// Node* n2=NULL;
				// n2=node_overlying(n2,vv2);n2=node_overlying(n2,vv1);n2=node_overlying(n2,v2);n2=node_overlying(n2,v1);
				re=node_overlying(re,vs);
			}
		}
		free_node(n1);
	}

	int_rb_tree_free_value(tree);int_rb_tree_free_value(tree1);
	int_rb_tree_free(tree);int_rb_tree_free(tree1);
	int_rb_tree_free(tree3);
	free_kdnode(kdtree);
	free(values);
	return re;
}

static void points_cloud_filled_tools_get_splited_boundaries(Node**n,Node* temp_n,template_v* nv)
{
	Node*ns[2]={NULL,NULL};
	int k=0;
	for(Node* nit=temp_n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v*v=(template_v*)(nit->value);
		if(v==nv)
		{
			k++;
		}
		ns[k%2]=node_pushback(ns[k%2],v);
	}
	ns[0]=node_reverse(ns[0]);ns[1]=node_reverse(ns[1]);
	*n=node_overlying(*n,ns[0]);
	if(ns[1]!=NULL)
	{
		*n=node_overlying(*n,ns[1]);
	}
}

static Node* points_cloud_filled_get_splited_boundaries_from_intersecting_boundaries(Node*boundaries,Node* intersections,Mesh* points,KD_Node*kdtree)
{
	if(boundaries==NULL)
	{
		return NULL;
	}
	Node* bound=node_copy(boundaries);
	Node* n2=NULL;
	for(Node* nit=intersections;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v**vs=(template_v**)(nit->value);
		double*l1=get_line_equation_from_two_points(vs[0]->point,vs[1]->point);
		double*l2=get_line_equation_from_two_points(vs[2]->point,vs[3]->point);
		double t=line_point_nearest_2_line(l1,l2);
		double p[3]={t*l1[0]+l1[3],t*l1[1]+l1[4],t*l1[2]+l1[5]};

		template_v* nv=points->create_vertexv(points,p,3);
		kd_tree_insert_vertex(kdtree,nv);	
		filled_normal_info_to_new_vertex(nv,0,kdtree);	

		Node* node=(Node*)(node_find(bound,vs[0]));
		node_insertback(node,nv);
		node=(Node*)(node_find(bound,vs[2]));
		node_insertback(node,nv);
		n2=node_overlying(n2,nv);
		free(l1);free(l2);
		free(vs);	
	}
	Node*re=node_overlying(NULL,bound);
	//Node* n1=node_overlying(NULL,bound);
	for(Node* nit=n2;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* nv=(template_v*)(nit->value);
		Node*n1=NULL;
		for(Node* nit1=re;nit1!=NULL;nit1=(Node*)(nit1->Next))
		{
			Node* temp_n=(Node*)(nit1->value);
			points_cloud_filled_tools_get_splited_boundaries(&n1,temp_n,nv);
			free_node(temp_n);
		}	
		free_node(re);
		re=n1;
	}
	return re;
}

// double * compute_plane_normal_of_jordan_curve1(double** v,int len )
// {
//     double* re=(double*)malloc(sizeof(double)*3); 
//     memset(re,0,sizeof(double)*3);
//     for(int i=0;i<len;i++)
//     {
//         int j=(i-1+len)%len,k=(i+1)%len;
//         double a[3]={v[i][0]-v[j][0],v[i][1]-v[j][1],v[i][2]-v[j][2]};
//         double b[3]={v[k][0]-v[i][0],v[k][1]-v[i][1],v[k][2]-v[i][2]};
//         if(normalize(a,3)==0)
//         {
//         	printf("zheli a cuowu\n");
//         }
//         if(normalize(b,3)==0)
//         {
//         	printf("zheli b cuowu\n");
//         }
        
//         double * n=out_product(a,b);
//         double temp=acos(inner_product(a,b,3));
//         if(normalize(n,3)==0)
//         {
//         	printf("zheli n cuowu\n");
//         }

//         re[0]+=temp*n[0];re[1]+=temp*n[1];re[2]+=temp*n[2];
//         //printf("n : %lf %lf %lf %lf\n",n[0],n[1],n[2],temp);
//         SAFE_FREE(n);
//     } 
//     return re;
// }

//double* compute_plane_normal_of_jordan_curve_from_node(Node*);
static inline double* points_cloud_filled_get_raw_normal_vector(Node*n)
{
	double* re=(double*)malloc(sizeof(double)*3);

	double* axis=compute_pca_result_using_eigen(n);
	re[0]=axis[6];re[1]=axis[7];re[2]=axis[8];
	free(axis);
	return re;
}
double*  points_cloud_filled_adjust_normal_vector(Node* n,double* raw_nor)
{
	
	double** values=(double**)malloc(sizeof(double*)*node_size(n));
	int i=0;
	template_v* v=(template_v*)(n->value);
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v1=(template_v*)(nit->value);		
		//double p[3]={v1->point[0],v1->point[1],v1->point[2] };

		double p[3]={v1->point[0]-v->point[0],v1->point[1]-v->point[1],v1->point[2]-v->point[2] };

		double temp=inner_product(p,raw_nor,3);
		p[0]=p[0]-temp*raw_nor[0];p[1]=p[1]-temp*raw_nor[1];p[2]=p[2]-temp*raw_nor[2];
		
		values[i]=(double*)malloc(sizeof(double)*3);
		//printf("id:%d vallus i:%lf %lf %lf\n ",v1->id,p[0],p[1],p[2]);
		values[i][0]=p[0];values[i][1]=p[1];values[i][2]=p[2];
		i++;
	}
	double *normal = compute_plane_normal_of_jordan_curve(values,i);

	for(int j=0;j<i;j++)
	{
		free(values[j]);
	}	
	free(values);
	if(normalize(normal,3)==0)
	{
		printf("cuowu\n");
	}
	printf(" normal : %lf %lf %lf %lf %lf %lf\n",raw_nor[0],raw_nor[1],raw_nor[2],normal[0],normal[1],normal[2]);
	return normal;	
} 
static inline double points_cloud_filled_tools_compute_angle(double*normal,template_v*v1,template_v* v2,template_v*v3)
{
	double p1[3]={v3->point[0]-v2->point[0],v3->point[1]-v2->point[1],v3->point[2]-v2->point[2]};
	double p2[3]={v1->point[0]-v2->point[0],v1->point[1]-v2->point[1],v1->point[2]-v2->point[2]};

	double temp1=inner_product(normal,p1,3);
	double temp2=inner_product(normal,p2,3);

	p1[0]=p1[0]-temp1*normal[0];p1[1]=p1[1]-temp1*normal[1];p1[2]=p1[2]-temp1*normal[2];	
	p2[0]=p2[0]-temp2*normal[0];p2[1]=p2[1]-temp2*normal[1];p2[2]=p2[2]-temp2*normal[2];
	
	normalize(p1,3);normalize(p2,3);		
	double* out_vector=out_product(p1,p2);
	double temp=inner_product(out_vector,normal,3);
	double temp3=inner_product(p1,p2,3);
	free(out_vector);
	return (temp<0?(2.0*M_PI-acos(temp3)):(acos(temp3)));

}


void points_cloud_filled_find_min_angle(Node* n,double* normal ,template_v**v1,template_v**v2,template_v**v3)
{
	double angle=0;
	template_v* v=(template_v*)(n->value);
	template_v* v_1=(template_v*)(node_reverse(n)->value);
	template_v* v_2=(template_v*)(((Node*)(n->Next))->value);
		
	angle=points_cloud_filled_tools_compute_angle(normal,v_1,v,v_2) ;	
	*v1=v_1;*v2=v;*v3=v_2;
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		v=(template_v*)(nit->value);
		v_1=(nit->Prev==NULL)?((template_v*)(node_reverse(nit)->value)):((template_v*)(((Node*)(nit->Prev))->value));
		v_2=(nit->Next==NULL)?((template_v*)(n->value)):((template_v*)(((Node*)(nit->Next))->value));
		
		double temp_angle= points_cloud_filled_tools_compute_angle(normal,v_1,v,v_2) ;
		//printf(" %d %d %d angle:%lf\n",v_1->id,v->id,v_2->id,temp_angle);
		if(temp_angle<angle)
		{
			angle=temp_angle;

			*v1=v_1;*v2=v;*v3=v_2;
		}
	}
}
static inline double* points_cloud_filled_tools_get_media_vector(double*normal,template_v* v1,template_v* v2)
{

	double p[3]={v2->point[0]-v1->point[0],v2->point[1]-v1->point[1],v2->point[2]-v1->point[2]};

	double temp=inner_product(normal,p,3);
	p[0]=p[0]-temp*normal[0];	p[1]=p[1]-temp*normal[1];	p[2]=p[2]-temp*normal[2];
	double * re=(double*)malloc(sizeof(double)*6);
	re[3]=v1->point[0]+p[0]/2.0;	re[4]=v1->point[1]+p[1]/2.0;	re[5]=v1->point[2]+p[2]/2.0;

	double * nor=out_product(p,normal);

	normalize(nor,3);
	re[0]=nor[0];re[1]=nor[1];re[2]=nor[2];

	free(nor);								
	return re;										

}
void points_cloud_filled_chuli_one_vertex(Mesh*points ,Node** n,double r,KD_Node*kdtree,double* normal1,template_v* v1,template_v*v2,template_v*v3)
{
	Normal_Boundary_Info*nbi=(Normal_Boundary_Info*)(v2->prop);
	double  normal[3]={nbi->normal[0],nbi->normal[1],nbi->normal[2]};
	if(inner_product(normal1,normal,3)<0)
	{
		normal[0]=-normal[0];normal[1]=-normal[1];normal[2]=-normal[2];
	}
	normalize(normal,3);
	normal[0]+=0.1*normal1[0];normal[1]+=0.1*normal1[1];normal[2]+=0.1*normal1[2];
	normalize(normal,3);

	double angle= points_cloud_filled_tools_compute_angle(normal,v1,v2,v3);	
	if(angle<M_PI/2.0)
	{
		printf("<M_PI/2.0\n");
		Node* node=node_find(*n,v2);

		*n=node_remove(*n,node);
		printf("end\n");
	}
	else if(angle<M_PI*135.0/180.0)
	{
		printf("<M_PI*135.0/180.0\n");
		Node* node=node_find(*n,v2);
		double p1[3]={v1->point[0]-v2->point[0], v1->point[1]-v2->point[1], v1->point[2]-v2->point[2] };
		double p2[3]={v3->point[0]-v2->point[0], v3->point[1]-v2->point[1], v3->point[2]-v2->point[2] };
		double p3[3]={p1[0]+p2[0], p1[1]+p2[1], p1[2]+p2[2]};
		double t=inner_product(p3,normal,3);
		p3[0]=p3[0]-t*normal[0]+v2->point[0];p3[1]=p3[1]-t*normal[1]+v2->point[1];p3[2]=p3[2]-t*normal[2]+v2->point[2];

		
		template_v*nv=points->create_vertexv(points,p3,3);
		printf("create %d\n",nv->id);
		node_insertback(node,nv);
		*n=node_remove(*n,node);
		kd_tree_insert_vertex(kdtree,nv);	
		filled_normal_info_to_new_vertex(nv,r,kdtree);
	}
	else if(angle<M_PI*200.0/180.0)
	{
		printf("<M_PI*200.0/180.0\n");
		double * l1=points_cloud_filled_tools_get_media_vector(normal,v2,v1);
		double * l2=points_cloud_filled_tools_get_media_vector(normal,v2,v3);
		Node*node=node_find(*n,v2);
		double t= distance_of_two_points(v2->point,v1->point,3)/2.0;
		double p[3]={t*l1[0]+l1[3],t*l1[1]+l1[4],t*l1[2]+l1[5]};
		template_v* nv1= points->create_vertexv(points,p,3) ;
		t= -distance_of_two_points(v2->point,v3->point,3)/2.0;
		p[0]=t*l2[0]+l2[3];p[1]=t*l2[1]+l2[4];p[2]=t*l2[2]+l2[5];
		template_v* nv2= points->create_vertexv(points,p,3);
		printf("create %d %d\n",nv1->id,nv2->id);

		kd_tree_insert_vertex(kdtree,nv1);kd_tree_insert_vertex(kdtree,nv2);
		filled_normal_info_to_new_vertex(nv1,r,kdtree);filled_normal_info_to_new_vertex(nv2,r,kdtree);
		if(node->Prev==NULL)
		{
			Node* node_n=node_reverse(*n);
			node_insertback(node_n,nv1);
		}
		else{
			Node* node_n=(Node*)(node->Prev);
			node_insertback(node_n,nv1);
		}	
		node_insertback(node,nv2);
		*n=node_remove(*n,node);
		free(l1);free(l2);	

	}
	else {
		printf(">M_PI*200.0/180.0\n");
		double * l1=points_cloud_filled_tools_get_media_vector(normal,v2,v1);
		double * l2=points_cloud_filled_tools_get_media_vector(normal,v2,v3);
		Node*node=node_find(*n,v2);
		double t= distance_of_two_points(v2->point,v1->point,3)/2.0;
		double p[3]={t*l1[0]+l1[3],t*l1[1]+l1[4],t*l1[2]+l1[5]};
		template_v* nv1= points->create_vertexv(points,p,3) ;
		t= -distance_of_two_points(v2->point,v3->point,3)/2.0;
		p[0]=t*l2[0]+l2[3];p[1]=t*l2[1]+l2[4];p[2]=t*l2[2]+l2[5];
		template_v* nv2= points->create_vertexv(points,p,3);
		printf("create %d %d\n",nv1->id,nv2->id);

		kd_tree_insert_vertex(kdtree,nv1);kd_tree_insert_vertex(kdtree,nv2);
		filled_normal_info_to_new_vertex(nv1,r,kdtree);filled_normal_info_to_new_vertex(nv2,r,kdtree);
		if(node->Prev==NULL)
		{
			Node* node_n=node_reverse(*n);
			node_insertback(node_n,nv1);
		}
		else{
			Node* node_n=(Node*)(node->Prev);
			node_insertback(node_n,nv1);
		}	
		node_insertback(node,nv2);
		//*n=node_remove(*n,node);
		free(l1);free(l2);	
	}
}

void points_cloud_filled_test(Mesh* points,Node**n,double r,KD_Node*kdtree)
{
	printf("begin points cloud filled  test\n");
	if(node_size(*n)<=3)
	{
		return;
	}
	double p1[3]={1.0,0.0,0.0},p2[3]={0.0,1.0,0.0};
	double* p3=out_product(p1,p2);
	printf("p3 :%lf %lf %lf\n",p3[0],p3[1],p3[2]);	
//**************************
	
	int i=0;
//**************************
	do{
//		break;
		printf("once\n");
		i++;
		points_cloud_filled_pre_handle(points,*n,r,kdtree);
		points_cloud_filled_pre_handle1(points,n,r,kdtree);	
		double* raw_nor= points_cloud_filled_get_raw_normal_vector(*n);

		printf("node size %d\n",node_size(*n));

		Node* n1=points_cloud_filled_intersections_of_boundary(*n,r,raw_nor);
		if(n1!=NULL)
		{
			Node* n2=points_cloud_filled_get_splited_boundaries_from_intersecting_boundaries(*n,n1,points,kdtree);
			for(Node* nit=n2;nit!=NULL;nit=(Node*)(nit->Next))
			{
				Node*temp_n=(Node*)(nit->value);
				if(node_size(temp_n)>3)
				{
					free_node(*n);
					*n=temp_n;
				}
				else{
					free_node(temp_n);
				}	
				//printf("%d \n",node_size(temp_n));
			}	
			printf("cuowu %d\n",node_size(n2));
			free_node(n2);	
			//break;
		}

		free_node(n1);	

		double* normal=points_cloud_filled_adjust_normal_vector(*n,raw_nor);
		free(raw_nor);	
		template_v* v1=NULL,*v2=NULL,*v3=NULL;	
		points_cloud_filled_find_min_angle(*n,normal,&v1,&v2,&v3);

		printf("v1 :%d v2: %d  v3:%d\n",v1->id,v2->id,v3->id);
		printf("v2: %lf %lf %lf\n",v2->point[0],v2->point[1],v2->point[2]);
		if(node_find(*n,v2)==NULL)
		{
			printf("not find v2\n");
		}
		points_cloud_filled_chuli_one_vertex(points,n,r,kdtree,normal,v1,v2,v3);
		free(normal);
		if(i>=2000)
		{
			break;
		}
	}while(node_size(*n)>3);
}





void test_show()
{

	Mesh mesh;
	Mesh_init(&mesh);
	//"cloud_result.off"
	//"rabbitCut.off"
	_ReadOff_(&mesh,"rabbitCut.off",3);

	mesh_adjust_mesh_to_unit_hexahedron(&mesh);	
	Node* n=NULL;
	for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
	{
		n=node_overlying(n,quote(vit));
	}
	double* box=mesh_compute_bounding_box_from_node(n);
	KD_Node* kdtree=create_kd_tree_from_mesh(&mesh);
    printf("box:%lf %lf %lf %lf %lf %lf\n",box[0],box[1],box[2],box[3],box[4],box[5]);
    
	Node* n1=points_cloud_compute_normal_and_get_boundary_points(n,kdtree,7);
	if(n1==NULL)
	{
		printf("is null\n");
	}
	double r=compute_visual_avarage_distance(n);
	printf("r :%lf\n",r);

	Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree1);
	Node* node_pos=compute_points_open_sets(n,tree1,r,kdtree);
	printf("node pos size:%d\n",node_size(node_pos));


	set_normal_vector_of_manifold_points(node_pos);

	Mesh* mesh1=test_connect_boundaries(n1,kdtree ,r );

	for(Node*nit=node_pos;nit!=NULL;nit=(Node*)(nit->Next))
	{
		points_open_set_free((Points_Open_Set*)(nit->value));	
	}	
	free_node(node_pos);
	for(auto it=tree1->begin(tree1);it.it!=NULL;it++)
	{
		free_node((Node*)(it.second));
	}
	int_rb_tree_free(tree1);


	Int_RB_Tree* tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree2);
	Node* node_edge= points_cloud_filled_get_one_boundaries(&mesh,mesh1,tree2);

	free_node(node_edge);
	node_edge= points_cloud_filled_get_one_boundaries(&mesh,mesh1,tree2);

	free_node(node_edge);
	node_edge= points_cloud_filled_get_one_boundaries(&mesh,mesh1,tree2);


	printf("num v :%d\n",mesh.num_v(&mesh));

	points_cloud_filled_test(&mesh,&node_edge,r,kdtree);
	printf("num v :%d\n",mesh.num_v(&mesh));

	int_rb_tree_free(tree2);

	Viewer_World_Manager vwm;
	viewer_world_manager_init(&vwm);
	Viewer_World *vw=vwm.create_world(&vwm,NULL);

	add_default_somethings(vw);

	//test_show_mesh_points(vw,&mesh);
	//*************************Points
	Node* n2=vw->create_something(vw,"Points");
 	Viewer_Something*vs=(Viewer_Something*)(n2->value);
 	Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
 	vp->Data_rows=mesh.num_v(&mesh);
 	vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
 	int i=0;
 	vp->color=(float*)malloc(sizeof(float)*4*vp->Data_rows);
 	memset(vp->color,1,sizeof(float)*4*vp->Data_rows);	
   	for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
   	{
   		for(int j=0;j<3;j++)
   		{
   			vp->Data[i*3+j]=quote(vit)->point[j];
   		}	
   		vp->color[i*4+0]=0.0;vp->color[i*4+1]=0.8;vp->color[i*4+2]=0.0;vp->color[i*4+3]=1.0;
   		Normal_Boundary_Info* nbi=(Normal_Boundary_Info*)(quote(vit)->prop);	
   		//if(nbi!=NULL&&nbi->is_boundary==1)
   		if(quote(vit)->id==151118)
   		{
   			printf("normal : %lf %lf %lf %lf %lf %lf\n",nbi->normal[0],nbi->normal[1],nbi->normal[2],quote(vit)->point[0],quote(vit)->point[1],quote(vit)->point[2]);

   			vp->color[i*4+0]=1.0;
   		}
   		i++;
   	} 
 	vp->pointsize=10.0;
	free_node(n2);
 	//******************************Edges

	n2=vw->create_something(vw,"Edges");

	vs=(Viewer_Something*)(n2->value);
	Viewer_Edges*ve=(Viewer_Edges*)(vs->evolution);
	ve->Data_rows=node_size(node_edge);
	ve->Data_index_rows=ve->Data_rows;
	ve->color_rows=ve->Data_index_rows;
	ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
	ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int )*2*ve->Data_index_rows);
	ve->color=(float*)malloc(sizeof(float)*4*ve->color_rows);
	ve->set_color(ve,1.0,0.0,0.0,1.0);
	i=0;
	for(Node* nit=node_edge;nit!=NULL;nit=(Node*)(nit->Next))
	{
		ve->Data_index[2*i+0]=i;ve->Data_index[2*i+1]=(i+1)%ve->Data_rows;
		for(int j=0;j<3;j++)
		{
			ve->Data[i*3+j]=((template_v*)(nit->value))->point[j];
		}	
		i++;
	}		

	free_node(n2);

	//**************************Edges
	n2=vw->create_something(vw,"Edges");
	vs=(Viewer_Something*)(n2->value);
	ve=(Viewer_Edges*)(vs->evolution);
	ve->Data_rows=2;
	ve->Data_index_rows=1;
	ve->color_rows=ve->Data_index_rows;
	ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
	ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);
	ve->color=(float*)malloc(sizeof(float)*4*ve->color_rows);
	ve->set_color(ve,1.0,0.0,0.0,1.0);
	ve->Data[0]=0.029886;ve->Data[1]=0.080902; ve->Data[2]=-0.541195;
	ve->Data[3]=ve->Data[0]-0.823743;ve->Data[4]=ve->Data[1]+0.412987; ve->Data[5]=ve->Data[2]-0.388445;
	ve->Data_index[0]=0;ve->Data_index[1]=1;

	free_node(n2);

	// n2=vw->create_something(vw,"Edges");
	// vs=(Viewer_Something*)(n2->value);
	// Viewer_Edges*ve=(Viewer_Edges*)(vs->evolution);
	// ve->Data_rows=mesh1->num_v(mesh1);
	// ve->Data_index_rows=mesh1->num_c(mesh1);
	// ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
	// ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);
	// ve->color_rows=ve->Data_index_rows;
	// ve->color= (float*)malloc(sizeof(float)*4* ve->color_rows);
	// ve->set_color(ve,1.0,0.0,0.0,1.0);			
	// i=0;
	// for(auto cit=mesh1->c_begin(mesh1);cit!=mesh1->c_end(mesh1);cit++)
	// {
	// 	ve->Data_index[i*2+0]=quote(cit)->vertices[0]->id;
	// 	ve->Data_index[i*2+1]=quote(cit)->vertices[1]->id;
		
	// 	for(int j=0;j<3;j++)
	// 	{
	// 		ve->Data[3*(quote(cit)->vertices[0]->id)+j]=quote(cit)->vertices[0]->point[j];
	// 		ve->Data[3*(quote(cit)->vertices[1]->id)+j]=quote(cit)->vertices[1]->point[j];	
	// 	}
	// 	i++;
	// }
	// free_node(n2);

	//****************************
	//************************Edges
	// n2=vw->create_something(vw,"Edges");
	// vs=(Viewer_Something*)(n2->value);
	// ve=(Viewer_Edges*)(vs->evolution);
	// ve->Data_rows=mesh.num_v(&mesh)*2;
	// ve->Data_index_rows=mesh.num_v(&mesh);

	// ve->color_rows=ve->Data_index_rows;
	// ve->Data=(float *)malloc(sizeof(float)*3*ve->Data_rows);	
	// ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);
	// ve->color=(float*)malloc(sizeof(float)*4*ve->color_rows);
	// ve->set_color(ve,0.0,0.0,0.0,1.0);	
	// i=0;
	// for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
	// {
	// 	for(int j=0;j<3;j++)
	// 	{
	// 		ve->Data[(i*2)*3+j]=quote(vit)->point[j];
	// 	}
	// 	Normal_Boundary_Info*nbi=(Normal_Boundary_Info*)(quote(vit)->prop);
	// 	for(int j=0;j<3;j++)
	// 	{
	// 		ve->Data[(i*2+1)*3+j]=quote(vit)->point[j]-nbi->normal[j]*0.1;
	// 	}	
	// 	ve->Data_index[i*2]=i*2;
	// 	ve->Data_index[i*2+1]=i*2+1;
	// 	//if(nbi->is_correct_normal==1)
	// 	{
	// 		ve->color[i*4+0]=1.0;
	// 	}
	// 	i++;
	// }	

	// free_node(n2);	
	//*********************

	Viewer_Opengl_Interpreter voi;
	viewer_opengl_interpreter_initn(&voi,&vwm);
	voi.interpreter(&voi);	

	free_node(n1); 
    printf("end box\n");
	free_kdnode(kdtree);
	free(box);
	free_node(n);
	viewer_opengl_interpreter_free(&voi);
	viewer_world_manager_free(&vwm);
	Mesh_free(&mesh);
}

void test_show1()
{
	Mesh mesh;
	Mesh_init(&mesh);
	//_ReadObj_(&mesh,"BAA.obj");
	_ReadOff_(&mesh,"cloud_result.off",3);
	//printf(" minafold require:%d m dimension:%d\n",mesh.manifold_require,mesh.dimension);	

	Mesh mesh1;
	Mesh_init(&mesh1);
	_ReadOff_(&mesh1,"quad.off",3);


	Viewer_World_Manager vwm;
	viewer_world_manager_init(&vwm);
	Viewer_World *vw=vwm.create_world(&vwm,NULL);

	add_default_somethings(vw);
	//test_show_mesh_cells(vw,&mesh);
	//test_show_mesh_cells(vw,&mesh1);
	Viewer_Something* vs=test_show_mesh_points(vw,&mesh);

	Interactive_Show_Info* isi=(Interactive_Show_Info*)malloc(sizeof(Interactive_Show_Info));
	interactive_show_info_init(isi);	
	Node* n=NULL;
	for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
	{
		n=node_overlying(n,quote(vit));
	}	
	isi->cur_n=n;
	isi->r=0.470978;
	isi->kdtree=create_kd_tree_from_mesh(&mesh);
	isi->m=&mesh;
	isi->vs=vs;
	isi->marked_vs= kd_tree_find_nearest_sphere((template_v*)(isi->cur_n->value),isi->r,isi->kdtree);	
	Node* n1=vw->create_something(vw,"Intera");
    vs=(Viewer_Something*)(n1->value);
    Viewer_Intera* vi=(Viewer_Intera*)(vs->evolution);

   	vi->key_callback=interactive_show_info_key_callback;
   	vi->representation=isi; 


	free_node(n1);
	//test_show_mesh_lines(vw,&mesh);

	Viewer_Opengl_Interpreter voi;
	viewer_opengl_interpreter_initn(&voi,&vwm);
	isi->voi=&voi;
	voi.interpreter(&voi);

	free_node(n);
	interactive_show_info_free(isi);

	viewer_opengl_interpreter_free(&voi);
	viewer_world_manager_free(&vwm);	
	Mesh_free(&mesh1);
	Mesh_free(&mesh);	

}

void compute_visual_distance()
{
	Mesh mesh;
	Mesh_init(&mesh);
	//"cloud_result.off"
	//"rabbitCut.off"
	_ReadOff_(&mesh,"cloud_result.off",3);

	Node* n=NULL;
	for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
	{
		n=node_overlying(n,quote(vit));
	}
	double* box=mesh_compute_bounding_box_from_node(n);
	KD_Node* kdtree=create_kd_tree_from_mesh(&mesh);

	printf("%lf %lf %lf \n",box[0],box[1],box[2]);

	printf("%lf %lf %lf \n",box[3],box[4],box[5]);	
	double min1=box[3]-box[0],min2=box[4]-box[1],min3=box[5]-box[2];
	double mid=min1>=min2?(min3>=min1?(min1):(min2>=min3?(min2):(min3))):(min3>=min2?(min2):(min1>=min3?(min1):(min3)));

	double vdis=mid/70.0;
	printf("vdis:%lf \n",vdis);
	double sum=0;
	for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_v* v=(template_v*)(nit->value);
		Node* n1=kd_tree_find_nearest_sphere(v,vdis,kdtree);
		sum+=node_size(n1);
		//printf("size :%d\n",node_size(n1));	
		free_node(n1);
	}	
	printf("avarage sum:%lf\n",sum/((double)(node_size(n))));

	free_kdnode(kdtree);

	free(box);
	free_node(n);	
	Mesh_free(&mesh);

}


int main()
{


	//compute_visual_distance();
	test_show();

	//test_show1();
	// Mesh mesh;
	// Mesh_init(&mesh);
	// mesh.dimension=1;mesh.manifold_require=1;mesh.simplex=1;
	
	// double point1[3]={0,0,0},point2[3]={1,0,0};
	// template_v* v1=mesh.create_vertexv(&mesh,point1,3);
	// template_v* v2=mesh.create_vertexv(&mesh,point2,3);
	// template_v* vs[2]={v2,v1};
	// mesh.create_cellv(&mesh,vs,2);
	// for(auto cit=mesh.c_begin(&mesh);cit!=mesh.c_end(&mesh);cit++)
	// {
	// 	printf("cid:%d\n",quote(cit)->id);
	// 	for(int i=0;i<quote(cit)->vertices_size;i++)
	// 	{
	// 		printf("%lf %lf %lf\n ",quote(cit)->vertices[i]->point[0],quote(cit)->vertices[i]->point[1],quote(cit)->vertices[i]->point[2]);
	// 	}
	// }
	// Mesh_free(&mesh);


	printf("end\n");
	return 0;
}
